基于 R 语言的科研论文绘图技巧详解(5)
点击下方公众号,回复资料分享,收获惊喜
简介
在查阅文献的过程中,看到了几幅非常不错的出版图,今天就跟着小编一起学习下,他们是怎么使用 R 绘制出来的。
今天主要介绍 第五幅图(E) —— 带置信区间的拟合曲线图及带误差项的散点图。这个图在科研绘图中较为常用,例如:绘制产品失效时间的散点图。
前四幅图的详细代码介绍可见:基于 R 语言的科研论文绘图技巧详解(4)、基于 R 语言的科研论文绘图技巧详解(3)基于 R 语言的科研论文绘图技巧详解(2)基于 R 语言的科研论文绘图技巧详解(1) 。最后一幅图会随后继续介绍,读者在学习过程中,可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列,更有利于读者学习和查阅。
那我们来看看,作者是怎么实现这个功能的吧,本文知识点较多,大家耐心学习,建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R[1] 中找到。
主要知识点
绘制散点图、丝带形状图; 绘制横向、纵向误差图; 学会小技巧:展示轴外部的图形。
绘图
加载包
首先加载一些需要使用到的包。
library(ggplot2) # Grammar of graphics
设置主题
接下来,为了方便起见,作者在绘图前设置好了主题,并将该函数命名为 my_theme
。
这一部分在第一篇推文给出,代码将在文末中完整代码给出。
手动修改大部分面板,具体可以参考本篇文章[2]。或者观看我在 B 站发布的《R 语言可视化教程》,里面也有一些简单主题设置介绍。
导入数据
首先使用 read.csv()
导入三个基因相关数据,最后合并到 data_E 数据框中。前 6 行数据如下所示。
data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
rep("gene b",nrow(data_Eb)),
rep("gene c",nrow(data_Ec))),
"t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
"C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
"pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
"neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
"pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
"neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)
head(data_E)
# gene t C pos_err neg_err pos_err_t neg_err_t
# 1 gene a 12 0.27 0.03 0.03 NA NA
# 2 gene a 24 0.19 0.02 0.02 NA NA
# 3 gene a 36 0.17 0.01 0.01 NA NA
# 4 gene a 48 0.13 0.04 0.04 NA NA
# 5 gene a 60 0.09 0.03 0.03 NA NA
# 6 gene a 72 0.10 0.01 0.01 NA NA
gene
是分类变量,t
表示时间将作为 x
轴。C
表示正常水平下的结果(小编对基因方面不是很懂,可能理解的不对)。pos_err
和 neg_err
是正向与反向误差。
自定义函数
首先,定义三个函数,以便稍后拟合散点图。这里的参数具体怎么得到的,作者没有介绍,小编不是很了解。或许是使用了某种参数估计的方法吧。
f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)
# 构造 f2 函数的数据
ribbon = data.frame(
"f2" = 0.3*exp(-t/60),
"t" = t
)
绘制图形
定义图形形状分别为实心正方形、圆形、三角形 manual_pch =c(15,16,17)
。之后简单绘制下散点图geom_point()
和 丝带形状的图 geom_ribbon()
。
注意:这里的
factor()
将 gene 变成离散形式。由于geom_ribbon()
中使用的数据不同,所以要重新指定下新数据。
manual_pch =c(15,16,17) 定义图形形状
ggplot(data=data_E) +
geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
fill="black",alpha=0.1)
在这个基础上使用 stat_function()
添加三个函数,两条虚线,一条实线。
stat_function(fun=f1, geom="line",linetype="dashed") +
stat_function(fun=f2, geom="line") +
stat_function(fun=f3, geom="line",linetype="dotted")
使用 geom_errorbar()
和 geom_errorbarh()
添加误差棒(纵向与横向)。注意内部参数不同:(ymin,ymax) 和 (xmin,xmax)。
geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), width=2) +
geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))
最后就是主题的设置了,包括运用了自定义的主题 my_theme()
。
scale_shape_manual(values=manual_pch) +
my_theme() + theme(legend.title=element_blank())+
theme(legend.position = c(0.9,0.95)) +
scale_x_continuous(expand = c(0, 0),
breaks = c(seq(0,96,12)),
limits = c(0,96)
) +
scale_y_continuous(expand = c(0, 0),
breaks = c(seq(0,0.4,0.05)),
limits = c(0,0.4)
) +
theme(
panel.grid.major = element_line("gray95", size = 0.1),
# putting label closer to axis bc exponent makes it bigger
axis.title.y = element_text(margin = margin(r = -9))
) +
xlab("time (h)") +
ylab(expression(paste("concentration (mmol",~cm^-3,")")))
但是,最右边几个点并没有显示清楚,这里作者运用了一个小技巧,将这些点清晰的呈现出来。来看看这个代码的效果:
coord_cartesian(clip = "off") # to allow for plotting outside axes
完整代码
# Panel E ----
library(ggplot2)
base_size = 12
my_theme <- function() {
theme(
aspect.ratio = 1,
axis.line =element_line(colour = "black"),
# shift axis text closer to axis bc ticks are facing inwards
axis.text.x = element_text(size = base_size*0.8, color = "black",
lineheight = 0.9,
margin=unit(c(0.3,0.3,0.3,0.3), "cm")),
axis.text.y = element_text(size = base_size*0.8, color = "black",
lineheight = 0.9,
margin=unit(c(0.3,0.3,0.3,0.3), "cm")),
axis.ticks = element_line(color = "black", size = 0.2),
axis.title.x = element_text(size = base_size,
color = "black",
margin = margin(t = -5)),
# t (top), r (right), b (bottom), l (left)
axis.title.y = element_text(size = base_size,
color = "black", angle = 90,
margin = margin(r = -5)),
axis.ticks.length = unit(-0.3, "lines"),
legend.background = element_rect(color = NA,
fill = NA),
legend.key = element_rect(color = "black",
fill = "white"),
legend.key.size = unit(0.5, "lines"),
legend.key.height =NULL,
legend.key.width = NULL,
legend.text = element_text(size = 0.6*base_size,
color = "black"),
legend.title = element_text(size = 0.6*base_size,
face = "bold",
hjust = 0,
color = "black"),
legend.text.align = NULL,
legend.title.align = NULL,
legend.direction = "vertical",
legend.box = NULL,
panel.background = element_rect(fill = "white",
color = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(size = base_size,
color = "black"),
)
}
data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
rep("gene b",nrow(data_Eb)),
rep("gene c",nrow(data_Ec))),
"t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
"C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
"pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
"neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
"pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
"neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)
head(data_E)
f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)
ribbon = data.frame(
"f2" = 0.3*exp(-t/60),
"t" = t
)
head(ribbon)
manual_pch =c(15,16,17) # available pch: type ?pch
panel_E <- ggplot(data=data_E) +
geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
fill="black",alpha=0.1) +
stat_function(fun=f1, geom="line",linetype="dashed") +
stat_function(fun=f2, geom="line") +
stat_function(fun=f3, geom="line",linetype="dotted") +
geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err),
width=2) +
geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))+
scale_shape_manual(values=manual_pch) +
my_theme() + theme(legend.title=element_blank())+
theme(legend.position = c(0.9,0.95)) +
scale_x_continuous(expand = c(0, 0),
breaks = c(seq(0,96,12)),
limits = c(0,96)
) +
scale_y_continuous(expand = c(0, 0),
breaks = c(seq(0,0.4,0.05)),
limits = c(0,0.4)
) +
theme(
panel.grid.major = element_line("gray95", size = 0.1),
# putting label closer to axis bc exponent makes it bigger
axis.title.y = element_text(margin = margin(r = -9))
) +
xlab("time (h)") +
ylab(expression(paste("concentration (mmol",~cm^-3,")"))) +
coord_cartesian(clip = "off") # to allow for plotting outside axes
panel_E
小编有话说
本文主要学到的知识点如下:
使用 geom_point()
绘制散点图,geom_ribbon()
绘制丝带形状图;使用 stat_function()
添加函数曲线;使用 geom_errorbar()
和geom_errorbarh
添加误差棒(纵向与横向);使用 coord_cartesian(clip = "off")
允许展示外轴的图形。
看完这篇文章,相信你的技能包又多了点东西。记得实操噢!看了不代表就会了~ 如果觉得内容有用的话,小编写的有心的话。给小编来杯咖啡吧!
参考资料
GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R: https://github.com/marco-meer/scifig_plot_examples_R
[2]文章: https://ggplot2.tidyverse.org/reference/theme.html
推荐: 可以保存以下照片,在b站扫该二维码,或者b站搜索【庄闪闪
】观看Rmarkdown系列的视频教程。Rmarkdown视频新增两节视频(写轮眼幻灯片制作)需要视频内的文档,可在公众号回复【rmarkdown
】
R沟通|Rmarkdown教程(4)
R沟通|Rmarkdown教程(3)
R沟通|Rmarkdown教程(2)